/* cddcore.c:  Core Procedures for cddlib
   written by Komei Fukuda, fukuda@ifor.math.ethz.ch
   Version 0.94, Aug. 4, 2005
*/

/* cddlib : C-library of the double description method for
   computing all vertices and extreme rays of the polyhedron 
   P= {x :  b - A x >= 0}.  
   Please read COPYING (GNU General Public Licence) and
   the manual cddlibman.tex for detail.
*/

using System;
using System.Collections.Generic;
using System.Text;

namespace CddSharp
{
    public partial class Cdd
    {

        public void CheckAdjacency(Cone cone,
            Ray RP1, Ray RP2, out bool adjacent)
        {
            Ray TempRay;
            bool localdebug = false;
            BitArray Face = null, Face1 = null;
            int last_m = 0;

            if (last_m != cone.m)
            {
                if (last_m > 0)
                {
                    Face = null; Face1 = null;
                }

                /*
                Face = new BitArray(cone.m);
                  Face1 = new BitArray(cone.m); 
                 * */

                last_m = cone.m;
            }

            if (debug) localdebug = true;

            adjacent = true;
            Face1 = RP1.ZeroSet.And(RP2.ZeroSet);
            Face = Face1.And(cone.AddedHalfspaces);

            if (Face.Count < cone.d - 2)
            {
                adjacent = false;
                if (localdebug)
                {
                    /*
                  fprintf(stderr,"non adjacent: set_card(face) %ld < %ld = cone.d.\n",
                    set_card(Face),cone.d);
                    */
                }
                return;
            }
            else if (cone.parent.NondegAssumed)
            {
                adjacent = true;
                return;
            }
            TempRay = cone.FirstRay;
            while (TempRay != null && adjacent)
            {
                if (TempRay != RP1 && TempRay != RP2)
                {
                    Face1 = TempRay.ZeroSet.And(cone.AddedHalfspaces);
                    if (Face.IsSubsetOf(Face1))
                        adjacent = false;
                }
                TempRay = TempRay.Next;
            }
        }

        public void Eliminate(Cone cone, Ray Ptr)
        {
            /*eliminate the record pointed by Ptr^.Next*/
            Ray TempPtr;

            TempPtr = Ptr.Next;
            Ptr.Next = Ptr.Next.Next;
            if (TempPtr == cone.FirstRay)   /*Update the first pointer*/
                cone.FirstRay = Ptr.Next;
            if (TempPtr == cone.LastRay)   /*Update the last pointer*/
                cone.LastRay = Ptr;

            /* Added, Marc Pfetsch 010219 */
            for (int j = 0; j < cone.d; j++)
                TempPtr.RRay[j] = 0.0;
            TempPtr.ARay = 0.0;

            TempPtr.RRay = null;          /* free the ray vector memory */
            TempPtr.ZeroSet = null;  /* free the ZeroSet memory */
            TempPtr = null;   /* free the Ray structure memory */
            cone.RayCount--;
        }

        public void SetInequalitySets(Cone cone)
        {

            cone.GroundSet.SetAll(false);
            cone.EqualitySet.SetAll(false);
            cone.NonequalitySet.SetAll(false);
            for (int i = 1; i <= (cone.parent.m); i++)
            {
                cone.GroundSet.Set(i, true);
                if (cone.parent.EqualityIndex[i] == 1) cone.EqualitySet.Set(i, true);
                if (cone.parent.EqualityIndex[i] == -1) cone.NonequalitySet.Set(i, true);
            }
        }


        public void AValue(out double val, int d_size, double[][] A, double[] p, int i)
        {
            /*return the ith component of the vector  A x p */
            int j;
            double x;

            x = 0.0;
            val = 0;
            /* Changed by Marc Pfetsch 010219 */

            for (j = 0; j < d_size; j++)
            {
                x = A[i - 1][j] * p[j];
                val += x;
            }
        }

        public void StoreRay1(Cone cone, double[] p, out bool feasible)
        {  /* Original ray storing routine when RelaxedEnumeration is false */
            int k, fii = cone.m + 1;
            double temp;
            Ray RR;
            bool localdebug = debug;

            temp = 0.0;
            RR = cone.LastRay;
            feasible = true;
            RR.ZeroSet = new BitArray(cone.m);
            for (int j = 0; j < cone.d; j++)
            {
                RR.RRay[j] = p[j];
            }
            for (int i = 1; i <= cone.m; i++)
            {
                k = cone.OrderVector[i];
                AValue(out temp, cone.d, cone.A, p, k);
                if (localdebug)
                {
                    /*
                  fprintf(stderr,"StoreRay1: AValue at row %ld =",k);
                  WriteNumber(stderr, temp);
                  fprintf(stderr,"\n");
                     * */
                }
                if (temp == 0.0)
                {
                    RR.ZeroSet.Set(k, true);
                    if (localdebug)
                    {
                        /*
                      fprintf(stderr,"recognized zero!\n");
                         * */
                    }
                }
                if (temp < 0)
                {
                    if (localdebug)
                    {
                        /*
                      fprintf(stderr,"recognized negative!\n");
                         * */
                    }
                    feasible = false;
                    if (fii > cone.m) fii = i;  /* the first violating inequality index */
                    if (localdebug)
                    {
                        /*
                      fprintf(stderr,"this ray is not feasible, neg comp = %ld\n", fii);
                      WriteNumber(stderr, temp);  fprintf(stderr,"\n");
                         * */
                    }
                }
            }
            RR.FirstInfeasIndex = fii;
            RR.feasible = feasible;
        }

        public void StoreRay2(Cone cone, double[] p,
            out bool feasible, out bool weaklyfeasible)
        /* Ray storing routine when RelaxedEnumeration is true.
            weaklyfeasible is true iff it is feasible with
            the strict_inequality conditions deleted. */
        {
            Ray RR;
            int k, fii = cone.m + 1;
            double temp;
            bool localdebug = debug;

            temp = 0.0;
            RR = cone.LastRay;
            if (debug) localdebug = true;
            feasible = true;
            weaklyfeasible = true;
            RR.ZeroSet = new BitArray(cone.m);
            for (int j = 0; j < cone.d; j++)
            {
                RR.RRay[j] = p[j];
            }
            for (int i = 1; i <= cone.m; i++)
            {
                k = cone.OrderVector[i];
                AValue(out temp, cone.d, cone.A, p, k);
                if (temp == 0.0)
                {
                    RR.ZeroSet.Set(k, true);
                    if (cone.parent.EqualityIndex[k] == -1)
                        feasible = false;  /* strict inequality required */
                }
                /*    if (temp < -zero){ */
                if (temp < 0)
                {
                    feasible = false;
                    if (fii > cone.m && cone.parent.EqualityIndex[k] >= 0)
                    {
                        fii = i;  /* the first violating inequality index */
                        weaklyfeasible = false;
                    }
                }
            }
            RR.FirstInfeasIndex = fii;
            RR.feasible = feasible;
        }

        public void AddRay(Cone cone, double[] p)
        {
            bool feasible, weaklyfeasible;

            if (cone.FirstRay == null)
            {
                cone.FirstRay = new Ray();
                cone.FirstRay.RRay = new double[cone.d];
                for (int j = 0; j < cone.d; j++)
                    cone.FirstRay.RRay[j] = 0.0;
                cone.FirstRay.ARay = 0.0;
                if (debug)
                {
                    // fprintf(stderr,"Create the first ray pointer\n");
                }
                cone.LastRay = cone.FirstRay;
                cone.ArtificialRay.Next = cone.FirstRay;
            }
            else
            {
                cone.LastRay.Next = new Ray();
                cone.LastRay.Next.RRay = new double[cone.d];
                for (int j = 0; j < cone.d; j++)
                    cone.LastRay.Next.RRay[j] = 0.0;
                cone.LastRay.Next.ARay = 0;
                if (debug)
                {
                    //   fprintf(stderr,"Create a new ray pointer\n");
                }
                cone.LastRay = cone.LastRay.Next;
            }
            cone.LastRay.Next = null;
            cone.RayCount++;
            cone.TotalRayCount++;
            if (debug)
            {
                /*
              if (cone.TotalRayCount % 100 == 0) {
                fprintf(stderr,"*Rays (Total, Currently Active, Feasible) =%8ld%8ld%8ld\n",
               cone.TotalRayCount, cone.RayCount, cone.FeasibleRayCount);
              }
                 * */
            }
            if (cone.parent.RelaxedEnumeration)
            {
                StoreRay2(cone, p, out feasible, out weaklyfeasible);
                if (weaklyfeasible) cone.WeaklyFeasibleRayCount++;
            }
            else
            {
                StoreRay1(cone, p, out feasible);
                if (feasible) cone.WeaklyFeasibleRayCount++;
                /* weaklyfeasible is equiv. to feasible in this case. */
            }
            if (!feasible) return;
            else
            {
                (cone.FeasibleRayCount)++;
            }
        }

        public void AddArtificialRay(Cone cone)
        {
            double[] zerovector;
            int d1;
            bool feasible;

            if (cone.d <= 0) d1 = 1; else d1 = cone.d;
            zerovector = new double[d1];
            if (cone.ArtificialRay != null)
            {
                /*
              fprintf(stderr,"Warning !!!  FirstRay in not nil.  Illegal Call\n");
              free(zerovector); /* 086 */
                //return;
            }
            cone.ArtificialRay = new Ray();
            cone.ArtificialRay.RRay = new double[d1];
            for (int j = 0; j < d1; j++)
                cone.ArtificialRay.RRay[j] = 0.0;
            cone.ArtificialRay.ARay = 0.0;

            if (debug)
            {
                /*
                fprintf(stderr,"Create the artificial ray pointer\n");
                 * */
            }

            cone.LastRay = cone.ArtificialRay;
            StoreRay1(cone, zerovector, out feasible);
            /* This stores a vector to the record pointed by cone.LastRay */
            cone.ArtificialRay.Next = null;
            for (int j = 0; j < d1; j++)
            {
                zerovector[j] = 0.0;
            }
            zerovector = null; /* 086 */
        }

        public void ConditionalAddEdge(Cone cone,
            Ray Ray1, Ray Ray2, Ray ValidFirstRay)
        {
            int it, it_row, fii1, fii2, fmin, fmax;
            bool adjacent, lastchance;
            Ray TempRay, Rmin, Rmax;
            Adjacency NewEdge;
            bool localdebug = false;
            BitArray ZSmin, ZSmax;
            BitArray Face, Face1;
            int last_m = 0;

            if (last_m != cone.m)
            {
                if (last_m > 0)
                {
                    Face = null; Face1 = null;
                }
                Face = new BitArray(cone.m);
                Face1 = new BitArray(cone.m);
                last_m = cone.m;
            }

            fii1 = Ray1.FirstInfeasIndex;
            fii2 = Ray2.FirstInfeasIndex;
            if (fii1 < fii2)
            {
                fmin = fii1; fmax = fii2;
                Rmin = Ray1;
                Rmax = Ray2;
            }
            else
            {
                fmin = fii2; fmax = fii1;
                Rmin = Ray2;
                Rmax = Ray1;
            }
            ZSmin = Rmin.ZeroSet;
            ZSmax = Rmax.ZeroSet;
            if (localdebug)
            {
                /*
              fprintf(stderr,"ConditionalAddEdge: FMIN = %ld (row%ld)   FMAX=%ld\n",
                fmin, cone.OrderVector[fmin], fmax);
                 * */
            }
            if (fmin == fmax)
            {
                /*
              if (localdebug) fprintf(stderr,"ConditionalAddEdge: equal FII value. No edge added\n");
                 * */
            }
            else if (ZSmax.Get(cone.OrderVector[fmin]))
            {
                /*
              if (localdebug) fprintf(stderr,"ConditionalAddEdge: No strong separation . No edge added\n");
                 * */
            }
            else
            {  /* the pair will be separated at the iteration fmin */
                lastchance = true;
                /* flag to check it will be the last chance to store the edge candidate */
                Face1 = ZSmax.And(ZSmin);
                cone.count_int++;
                if (localdebug)
                {
                    /*
                  fprintf(stderr,"Face: ");
                  for (it=1; it<=cone.m; it++) {
                    it_row=cone.OrderVector[it];
                    if (set_member(it_row, Face1)) fprintf(stderr,"%ld ",it_row);
                  }
                  fprintf(stderr,"\n");
                     * */
                }
                for (it = cone.Iteration + 1; it < fmin && lastchance; it++)
                {
                    it_row = cone.OrderVector[it];
                    if (cone.parent.EqualityIndex[it_row] >= 0 && Face1.Get(it_row))
                    {
                        lastchance = false;
                        cone.count_int_bad++;
                        if (localdebug)
                        {
                            /*
                          fprintf(stderr,"There will be another chance iteration %ld (row %ld) to store the pair\n", it, it_row);
                             * */
                        }
                    }
                }
                if (lastchance)
                {
                    adjacent = true;
                    (cone.count_int_good)++;
                    /* adjacent checking */
                    Face = Face1.And(cone.AddedHalfspaces);
                    if (localdebug)
                    {
                        /*
                      fprintf(stderr,"Check adjacency\n");
                      fprintf(stderr,"AddedHalfspaces: "); set_fwrite(stderr,cone.AddedHalfspaces);
                      fprintf(stderr,"Face: ");
                      for (it=1; it<=cone.m; it++) {
                        it_row=cone.OrderVector[it];
                        if (set_member(it_row, Face)) fprintf(stderr,"%ld ",it_row);
                      }
                      fprintf(stderr,"\n");
                         * */
                    }
                    if (Face.Count < cone.d - 2)
                    {
                        adjacent = false;
                    }
                    else if (cone.parent.NondegAssumed)
                    {
                        adjacent = true;
                    }
                    else
                    {
                        TempRay = ValidFirstRay;  /* the first ray for adjacency checking */
                        while (TempRay != null && adjacent)
                        {
                            if (TempRay != Ray1 && TempRay != Ray2)
                            {
                                Face1 = TempRay.ZeroSet.And(cone.AddedHalfspaces);
                                if (Face.IsSubsetOf(Face1))  //*****
                                {
                                    if (localdebug)
                                    {
                                        /*
                                        set_fwrite(stderr,Face1);
                                         * */
                                    }
                                    adjacent = false;
                                }
                            }
                            TempRay = TempRay.Next;
                        }
                    }
                    if (adjacent)
                    {
                        if (localdebug)
                        {
                            /*
                            fprintf(stderr,"The pair is adjacent and the pair must be stored for iteration %ld (row%ld)\n",
                          fmin, cone.OrderVector[fmin]);
                             * */
                        }
                        NewEdge = new Adjacency();
                        NewEdge.Ray1 = Rmax;  /* save the one remains in iteration fmin in the first */
                        NewEdge.Ray2 = Rmin;  /* save the one deleted in iteration fmin in the second */
                        NewEdge.Next = null;
                        (cone.EdgeCount)++;
                        (cone.TotalEdgeCount)++;
                        if (cone.Edges[fmin] == null)
                        {
                            cone.Edges[fmin] = NewEdge;
                            if (localdebug)
                            {
                                /*
                                fprintf(stderr,"Create a new edge list of %ld\n", fmin);
                                 * */
                            }
                        }
                        else
                        {
                            NewEdge.Next = cone.Edges[fmin];
                            cone.Edges[fmin] = NewEdge;
                        }
                    }
                }
            }
        }

        public void CreateInitialEdges(Cone cone)
        {
            Ray Ptr1, Ptr2;
            int fii1, fii2;
            int count = 0;
            bool adj, localdebug = false;

            cone.Iteration = cone.d;  /* CHECK */
            if (cone.FirstRay == null || cone.LastRay == null)
            {
                /* fprintf(stderr,"Warning:  CreateInitialEdges called with null pointer(s)\n"); */
                return;
            }
            Ptr1 = cone.FirstRay;
            while (Ptr1 != cone.LastRay && Ptr1 != null)
            {
                fii1 = Ptr1.FirstInfeasIndex;
                Ptr2 = Ptr1.Next;
                while (Ptr2 != null)
                {
                    fii2 = Ptr2.FirstInfeasIndex;
                    count++;
                    if (localdebug)
                    {
                        /*
                        fprintf(stderr, " CreateInitialEdges: edge %ld \n", count);
                         * */
                    }
                    CheckAdjacency(cone, Ptr1, Ptr2, out adj);
                    if (fii1 != fii2 && adj)
                        ConditionalAddEdge(cone, Ptr1, Ptr2, cone.FirstRay);
                    Ptr2 = Ptr2.Next;
                }
                Ptr1 = Ptr1.Next;
            }

        }


        public void UpdateEdges(Cone cone, Ray RRbegin, Ray RRend)
        /* This procedure must be called after the ray list is sorted
           by EvaluateARay2 so that FirstInfeasIndex's are monotonically
           increasing.
        */
        {
            Ray Ptr1, Ptr2begin, Ptr2;
            int fii1;
            bool ptr2found, quit, localdebug = false;
            int count = 0, pos1, pos2;
            double workleft, prevworkleft = 110.0, totalpairs;

            totalpairs = (cone.ZeroRayCount - 1.0) * (cone.ZeroRayCount - 2.0) + 1.0;
            Ptr2begin = null;
            if (RRbegin == null || RRend == null)
            {
                if (true)
                {
                    /*
                    fprintf(stderr,"Warning: UpdateEdges called with null pointer(s)\n");
                     * */
                }
                return;
            }
            Ptr1 = RRbegin;
            pos1 = 1;
            do
            {
                ptr2found = false;
                quit = false;
                fii1 = Ptr1.FirstInfeasIndex;
                pos2 = 2;
                for (Ptr2 = Ptr1.Next; !ptr2found && !quit; Ptr2 = Ptr2.Next, pos2++)
                {
                    if (Ptr2.FirstInfeasIndex > fii1)
                    {
                        Ptr2begin = Ptr2;
                        ptr2found = true;
                    }
                    else if (Ptr2 == RRend) quit = true;
                }
                if (ptr2found)
                {
                    quit = false;
                    for (Ptr2 = Ptr2begin; !quit; Ptr2 = Ptr2.Next)
                    {
                        count++;
                        if (localdebug)
                        {
                            /*
                            fprintf(stderr,"UpdateEdges: edge %ld \n",count);
                             * */
                        }
                        ConditionalAddEdge(cone, Ptr1, Ptr2, RRbegin);
                        if (Ptr2 == RRend || Ptr2.Next == null) quit = true;
                    }
                }
                Ptr1 = Ptr1.Next;
                pos1++;
                workleft = 100.0 * (cone.ZeroRayCount - pos1) * (cone.ZeroRayCount - pos1 - 1.0) / totalpairs;
                if (cone.ZeroRayCount >= 500 && debug && pos1 % 10 == 0 && prevworkleft - workleft >= 10)
                {
                    /*
                  fprintf(stderr,"*Work of iteration %5ld(/%ld): %4ld/%4ld => %4.1f%% left\n",
                     cone.Iteration, cone.m, pos1, cone.ZeroRayCount, workleft);
                     * */
                    prevworkleft = workleft;
                }
            } while (Ptr1 != RRend && Ptr1 != null);

        }

#if foobar
void FreeDDMemory0(Cone cone)
{

  Ray Ptr, PrevPtr;
  int count;
  int j;
  bool localdebug=false;
  
  /* THIS SHOULD BE REWRITTEN carefully */
  PrevPtr=cone.ArtificialRay;
  if (PrevPtr!=null){
    count=0;
    for (Ptr=cone.ArtificialRay.Next; Ptr!=null; Ptr=Ptr.Next){
      /* Added Marc Pfetsch 2/19/01 */
      for (j=0;j < cone.d;j++)
clear(PrevPtr.Ray[j]);
      clear(PrevPtr.ARay);

      free(PrevPtr.Ray);
      free(PrevPtr.ZeroSet);
      free(PrevPtr);
      count++;
      PrevPtr=Ptr;
    };
    cone.FirstRay=null;
    /* Added Marc Pfetsch 010219 */
    for (j=0;j < cone.d;j++)
clear(cone.LastRay.Ray[j]);
    clear(cone.LastRay.ARay);

    free(cone.LastRay.Ray);
    cone.LastRay.Ray = null;
    set_free(cone.LastRay.ZeroSet);
    cone.LastRay.ZeroSet = null;
    free(cone.LastRay);
    cone.LastRay = null;
    cone.ArtificialRay=null;
    if (localdebug) fprintf(stderr,"%ld ray storage spaces freed\n",count);
  }
/* must add (by Sato) */
  free(cone.Edges);
  
  set_free(cone.GroundSet); 
  set_free(cone.EqualitySet); 
  set_free(cone.NonequalitySet); 
  set_free(cone.AddedHalfspaces); 
  set_free(cone.WeaklyAddedHalfspaces); 
  set_free(cone.InitialHalfspaces);
  free(cone.InitialRayIndex);
  free(cone.OrderVector);
  free(cone.newcol);

/* Fixed by Shawn Rusaw.  Originally it was cone.d instead of cone.d_alloc */
  FreeBmatrix(cone.d_alloc,cone.B);
  FreeBmatrix(cone.d_alloc,cone.Bsave);

/* Fixed by Marc Pfetsch 010219*/
  FreeAmatrix(cone.m_alloc,cone.d_alloc,cone.A);
  cone.A = null;

  free(cone);


}

void FreeDDMemory(PolyhedraPtr poly)
{
  FreeDDMemory0(poly.child);
  poly.child=null;
}

void FreePolyhedra(PolyhedraPtr poly)
{
  int i;

  if ((poly).child != null) FreeDDMemory(poly);
  FreeAmatrix((poly).m_alloc,poly.d_alloc, poly.A);
  FreeArow((poly).d_alloc,(poly).c);
  free((poly).EqualityIndex);
  if (poly.AincGenerated){
    for (i=1; i<=poly.m1; i++){
      set_free(poly.Ainc[i-1]);
    }
    free(poly.Ainc);
    set_free(poly.Ared);
    set_free(poly.Adom);
    poly.Ainc=null;
  }

  free(poly);
}
#endif

        public void Normalize(int d_size, double[] V)
        {
            int jmin = 0;
            double temp, min;
            bool nonzerofound = false;

            if (d_size > 0)
            {
                min = 0.0;
                temp = 0.0;
                min = Math.Abs(V[0]);
                jmin = 0; /* set the minmizer to 0 */
                if (min > 0) nonzerofound = true;
                for (int j = 1; j < d_size; j++)
                {
                    temp = Math.Abs(V[j]);
                    if (temp > 0)
                    {
                        if (!nonzerofound || temp < min)
                        {
                            nonzerofound = true;
                            min = temp;
                            jmin = j;
                        }
                    }
                }
                if (min > 0)
                {
                    for (int j = 0; j < d_size; j++)
                        V[j] = V[j] / min;
                }
                //clear(min); clear(temp);
            }
        }


        public void ZeroIndexSet(int m_size, int d_size, double[][] A, double[] x, BitArray ZS)
        {
            int i;
            double temp;

            /* Changed by Marc Pfetsch 010219 */
            temp = 0;
            ZS.SetAll(false);
            for (i = 1; i <= m_size; i++)
            {
                AValue(out temp, d_size, A, x, i);
                if (temp == 0.0) ZS.Set(i, true);
            }

            /* Changed by Marc Pfetsch 010219 */
            //clear(temp);
        }

        public void CopyBmatrix(int d_size, double[][] T, double[][] TCOPY)
        {


            for (int i = 0; i < d_size; i++)
            {
                for (int j = 0; j < d_size; j++)
                {
                    TCOPY[i][j] = T[i][j];
                }
            }
        }


        public void CopyArow(double[] acopy, double[] a, int d)
        {
            for (int j = 0; j < d; j++)
            {
                acopy[j] = a[j];
            }
        }

        public void CopyNormalizedArow(double[] acopy, double[] a, int d)
        {
            CopyArow(acopy, a, d);
            Normalize(d, acopy);
        }

        public void CopyAmatrix(double[][] Acopy, double[][] A, int m, int d)
        {
            int i;

            for (i = 0; i < m; i++)
            {
                CopyArow(Acopy[i], A[i], d);
            }
        }

        public void CopyNormalizedAmatrix(double[][] Acopy, double[][] A, int m, int d)
        {
            int i;

            for (i = 0; i < m; i++)
            {
                CopyNormalizedArow(Acopy[i], A[i], d);
            }
        }

        public void PermuteCopyAmatrix(double[][] Acopy, double[][] A, int m, int d, int[] roworder)
        {
            int i;

            for (i = 1; i <= m; i++)
            {
                CopyArow(Acopy[i - 1], A[roworder[i] - 1], d);
            }
        }

        public void PermutePartialCopyAmatrix(double[][] Acopy, double[][] A, int m, int d, int[] roworder, int p, int q)
        {
            /* copy the rows of A whose roworder is positive.  roworder[i] is the row index of the copied row. */
            int i, k;

            k = 0;
            for (i = 1; i <= m; i++)
            {
                if (roworder[i] > 0) CopyArow(Acopy[roworder[i] - 1], A[i - 1], d);
            }
        }


        public void InitializeArow(int d, out double[] a)
        {
            /*
          int j;

          if (d>0) *a=(double*) calloc(d,sizeof(double));
          for (j = 0; j < d; j++) {
              init((*a)[j]);
          }
             * */
            a = new double[d];
        }

        public void InitializeAmatrix(int m, int d, out double[][] A)
        {

            A = new double[m][];
            for (int i = 0; i < m; i++)
            {
                InitializeArow(d, out A[i]);
            }
        }

        /*
void FreeAmatrix(int m,int d,double[][] A)
{
  int i;
  int j;

  for (i = 0; i < m; i++) {
    for (j = 0; j < d; j++) {
      clear(A[i][j]);
    }
  }
  if (A!=null) {
    for (i = 0; i < m; i++) {
      free(A[i]);
    }
    free(A);
  }
}

void FreeArow(int d, double[] a)
{
  int j;

  for (j = 0; j < d; j++) {
    clear(a[j]);
  }
  free(a);
}
*/

        public void InitializeBmatrix(int d, out double[][] B)
        {


            B = new double[d][];

            for (int j = 0; j < d; j++)
            {
                B[j] = new double[d];
            }
            /*for (i = 0; i < d; i++) {
              for (j = 0; j < d; j++) {
                init((*B)[i][j]);
              }
            }*/
        }

        /*
void FreeBmatrix(int d,double[][] B)
{
  int i,j;

  for (i = 0; i < d; i++) {
    for (j = 0; j < d; j++) {
      clear(B[i][j]);
    }
  }
  if (B!=null) {
    for (j = 0; j < d; j++) {
      free(B[j]);
    }
    free(B);
  }
}
        */


        public SetFamily CreateSetFamily(int fsize, int ssize)
        {
            SetFamily F;
            int i, f0, f1, s0, s1;

            if (fsize <= 0)
            {
                f0 = 0; f1 = 1;
                /* if fsize<=0, the fsize is set to zero and the created size is one */
            }
            else
            {
                f0 = fsize; f1 = fsize;
            }
            if (ssize <= 0)
            {
                s0 = 0; s1 = 1;
                /* if ssize<=0, the ssize is set to zero and the created size is one */
            }
            else
            {
                s0 = ssize; s1 = ssize;
            }

            F = new SetFamily();
            F.set = new BitArray[f1];
            for (i = 0; i < f1; i++)
            {
                F.set[i] = new BitArray(s1);
            }
            F.famsize = f0;
            F.setsize = s0;
            return F;
        }

        /*
void FreeSetFamily(SetFamilyPtr F)
{
  int i,f1;

  if (F!=null){
    if (F.famsize<=0) f1=1; else f1=F.famsize; 
      //* the smallest created size is one 
    for (i=0; i<f1; i++) {
      set_free(F.set[i]);
    }
    free(F.set);
    free(F);
  }
}
         * */

        public Matrix CreateMatrix(int m_size, int d_size)
        {
            Matrix M;
            int m0, m1;
            int d0, d1;

            if (m_size <= 0)
            {
                m0 = 0; m1 = 1;
                /* if m_size <=0, the number of rows is set to zero, the actual size is 1 */
            }
            else
            {
                m0 = m_size; m1 = m_size;
            }
            if (d_size <= 0)
            {
                d0 = 0; d1 = 1;
                /* if d_size <=0, the number of cols is set to zero, the actual size is 1 */
            }
            else
            {
                d0 = d_size; d1 = d_size;
            }
            M = new Matrix();
            InitializeAmatrix(m1, d1, out M.matrix);
            InitializeArow(d1, out M.rowvec);
            M.rowsize = m0;
            M.linset = new BitArray(m1);
            M.colsize = d0;
            M.objective = LPObjectiveType.LPnone;
            M.numbtype = NumberType.Unknown;
            M.representation = RepresentationType.Unspecified;
            return M;
        }

        /*
void FreeMatrix(Matrix M)
{
  int m1;
  int d1;

  if (M!=null) {
    if (M.rowsize<=0) m1=1; else m1=M.rowsize;
    if (M.colsize<=0) d1=1; else d1=M.colsize;
    if (M!=null) {
      FreeAmatrix(m1,d1,M.matrix);
      FreeArow(d1,M.rowvec);
      set_free(M.linset);
      free(M);
    }
  }
}
         * */

        public void SetToIdentity(int d_size, double[][] T)
        {
            for (int j1 = 1; j1 <= d_size; j1++)
            {
                for (int j2 = 1; j2 <= d_size; j2++)
                {
                    if (j1 == j2)
                        T[j1 - 1][j2 - 1] = 1.0;
                    else
                        T[j1 - 1][j2 - 1] = 0.0;
                }
            }
        }

        public void ColumnReduce(Cone cone)
        {
            int j, j1 = 0;
            int i;
            bool localdebug = false;

            for (j = 1; j <= cone.d; j++)
            {
                if (cone.InitialRayIndex[j] > 0)
                {
                    j1 = j1 + 1;
                    if (j1 < j)
                    {
                        for (i = 1; i <= cone.m; i++) cone.A[i - 1][j1 - 1] = cone.A[i - 1][j - 1];
                        cone.newcol[j] = j1;
                        if (localdebug)
                        {
                            /*
                          fprintf(stderr,"shifting the column %ld to column %ld\n", j, j1);
                             * */
                        }
                        /* shifting forward */
                    }
                }
                else
                {
                    cone.newcol[j] = 0;
                    if (localdebug)
                    {
                        /*
                      fprintf(stderr,"a generator (or an equation) of the linearity space: ");
                      for (i=1; i<=cone.d; i++) WriteNumber(stderr, cone.B[i-1][j-1]);
                      fprintf(stderr,"\n");
                         * */
                    }
                }
            }
            cone.d = j1;  /* update the dimension. cone.d_orig remembers the old. */
            CopyBmatrix(cone.d_orig, cone.B, cone.Bsave);
            /* save the dual basis inverse as Bsave.  This matrix contains the linearity space generators. */
            cone.ColReduced = true;
        }

        //rowbasis and colbasis might be out parameters?
        public int MatrixRank(Matrix M, BitArray ignoredrows, BitArray ignoredcols, out BitArray rowbasis, out BitArray colbasis)
        {
            bool stop, chosen, localdebug = debug;
            BitArray NopivotRow, PriorityRow;
            BitArray ColSelected;
            double[][] B;
            int[] roworder;
            int r;
            int s;
            int rank;

            rank = 0;
            stop = false;
            ColSelected = new BitArray(M.colsize);
            NopivotRow = new BitArray(M.rowsize);
            rowbasis = new BitArray(M.rowsize);
            colbasis = new BitArray(M.colsize);
            PriorityRow = new BitArray(M.rowsize);
            ignoredrows.CopyTo(NopivotRow);
            ignoredcols.CopyTo(ColSelected);
            InitializeBmatrix(M.colsize, out B);
            SetToIdentity(M.colsize, B);
            roworder = new int[M.rowsize + 1];
            for (int j = 0; j < M.rowsize; j++)
            {
                roworder[j + 1] = j + 1;
            }

            do
            {   /* Find a set of rows for a basis */
                SelectPivot2(M.rowsize, M.colsize, M.matrix, B, RowOrderType.MinIndex, roworder,
                 PriorityRow, M.rowsize, NopivotRow, ColSelected, out r, out s, out chosen);
                if (debug && chosen)
                {
                    /*
                  fprintf(stderr,"Procedure MatrixRank: pivot on (r,s) =(%ld, %ld).\n", r, s);
                     * */
                }
                if (chosen)
                {
                    NopivotRow.Set(r, true);
                    rowbasis.Set(r, true);
                    ColSelected.Set(s, true);
                    colbasis.Set(s, true);
                    rank++;
                    GaussianColumnPivot(M.rowsize, M.colsize, M.matrix, B, r, s);
                    if (localdebug)
                    {
                        /*
                        WriteBmatrix(stderr,M.colsize,B);
                         * */
                    }
                }
                else
                {
                    stop = true;
                }
                if (rank == M.colsize) stop = true;
            } while (!stop);
            /*
          FreeBmatrix(M.colsize,B);
          free(roworder);
          set_free(ColSelected);
          set_free(NopivotRow);
          set_free(PriorityRow);
             * */
            return rank;
        }


        public void FindBasis(Cone cone, out int rank)
        {
            bool stop, chosen, localdebug = debug;
            BitArray NopivotRow;
            BitArray ColSelected;
            int r;
            int j, s;

            rank = 0;
            stop = false;
            for (j = 0; j <= cone.d; j++) cone.InitialRayIndex[j] = 0;
            cone.InitialHalfspaces.SetAll(false); //set_emptyset
            ColSelected = new BitArray(cone.d);
            NopivotRow = new BitArray(cone.m);
            cone.NonequalitySet.CopyTo(NopivotRow);
            SetToIdentity(cone.d, cone.B);
            do
            {   /* Find a set of rows for a basis */
                SelectPivot2(cone.m, cone.d, cone.A, cone.B, cone.HalfspaceOrder, cone.OrderVector,
                 cone.EqualitySet, cone.m, NopivotRow, ColSelected, out r, out s, out chosen);
                if (debug && chosen)
                {
                    /*
                  fprintf(stderr,"Procedure FindBasis: pivot on (r,s) =(%ld, %ld).\n", r, s);
                     * */
                }
                if (chosen)
                {
                    cone.InitialHalfspaces.Set(r, true);
                    NopivotRow.Set(r, true);
                    ColSelected.Set(s, true);
                    cone.InitialRayIndex[s] = r;    /* cone.InitialRayIndex[s] stores the corr. row index */
                    rank++;
                    GaussianColumnPivot(cone.m, cone.d, cone.A, cone.B, r, s);
                    if (localdebug)
                    {
                        /*
                        WriteBmatrix(stderr,cone.d,cone.B);
                         * */
                    }
                }
                else
                {
                    stop = true;
                }
                if (rank == cone.d) stop = true;
            } while (!stop);

            //set_free(ColSelected);
            //set_free(NopivotRow);
        }


        public void FindInitialRays(Cone cone, out bool found)
        {
            BitArray CandidateRows;
            int i;
            int rank;
            RowOrderType roworder_save = RowOrderType.LexMin;

            found = false;
            CandidateRows = new BitArray(cone.m);
            if (cone.parent.InitBasisAtBottom == true)
            {
                roworder_save = cone.HalfspaceOrder;
                cone.HalfspaceOrder = RowOrderType.MaxIndex;
                cone.PreOrderedRun = false;
            }
            else cone.PreOrderedRun = true;
            if (debug)
            {
                /*
                WriteBmatrix(stderr, cone.d, cone.B);
                 * */
            }
            for (i = 1; i <= cone.m; i++)
                if (!cone.NonequalitySet.Get(i)) CandidateRows.Set(i, true);
            /*all rows not in NonequalitySet are candidates for initial cone*/
            FindBasis(cone, out rank);
            if (debug)
            {
                /*
                WriteBmatrix(stderr, cone.d, cone.B);
            fprintf(stderr,"FindInitialRays: rank of double[][] = %ld\n", rank);
                 * */
            }
            cone.LinearityDim = cone.d - rank;
            if (debug)
            {
                /*
                fprintf(stderr,"Linearity Dimension = %ld\n", cone.LinearityDim);
                 * */
            }
            if (cone.LinearityDim > 0)
            {
                ColumnReduce(cone);
                FindBasis(cone, out rank);
            }
            //if (!set_subset(cone.EqualitySet,cone.InitialHalfspaces)) {
            if (cone.EqualitySet.IsSubsetOf(cone.InitialHalfspaces))     ///***** subset
            {
                if (debug)
                {
                    /*
                  fprintf(stderr,"Equality set is dependent. Equality Set and an initial basis:\n");
                  set_fwrite(stderr,cone.EqualitySet);
                  set_fwrite(stderr,cone.InitialHalfspaces);
                     * */
                };
            }
            found = true;
            //set_free(CandidateRows);
            if (cone.parent.InitBasisAtBottom == true)
            {
                cone.HalfspaceOrder = roworder_save;
            }
            if (cone.HalfspaceOrder == RowOrderType.MaxCutoff ||
                cone.HalfspaceOrder == RowOrderType.MinCutoff ||
                cone.HalfspaceOrder == RowOrderType.MixCutoff)
            {
                cone.PreOrderedRun = false;
            }
            else cone.PreOrderedRun = true;
        }

        public void CheckEquality(int d_size, Ray RP1, Ray RP2, out bool equal)
        {
            int j;

            if (debug)
            {
                /*
              fprintf(stderr,"Check equality of two rays\n");
                 * */
            }
            equal = true;
            j = 1;
            while (j <= d_size && equal)
            {
                if (RP1.RRay[j - 1] != RP2.RRay[j - 1])
                    equal = false;
                j++;
            }
            if (equal)
            {
                /*
              fprintf(stderr,"Equal records found !!!!\n");
                 * */
            }
        }

        public void CreateNewRay(Cone cone,
            Ray Ptr1, Ray Ptr2, int ii)
        {
            /*Create a new ray by taking a linear combination of two rays*/
            int j;
            double a1, a2, v1, v2;
            double[] NewRay = null;
            int last_d = 0;
            bool localdebug = debug;

            a1 = 0.0; a2 = 0.0; v1 = 0.0; v2 = 0.0;
            if (last_d != cone.d)
            {
                if (last_d > 0)
                {
                    //for (j=0; j<last_d; j++) clear(NewRay[j]);
                    //free(NewRay);
                }
                NewRay = new double[cone.d];
                for (j = 0; j < cone.d; j++) NewRay[j] = 0.0;
                last_d = cone.d;
            }

            AValue(out a1, cone.d, cone.A, Ptr1.RRay, ii);
            AValue(out a2, cone.d, cone.A, Ptr2.RRay, ii);
            if (localdebug)
            {
                /*
              fprintf(stderr,"CreatNewRay: Ray1 ="); WriteArow(stderr, Ptr1.Ray, cone.d);
              fprintf(stderr,"CreatNewRay: Ray2 ="); WriteArow(stderr, Ptr2.Ray, cone.d);
                 * */
            }
            v1 = Math.Abs(a1);
            v2 = Math.Abs(a2);
            if (localdebug)
            {
                /*
              fprintf(stderr,"AValue1 and ABS");  WriteNumber(stderr,a1); WriteNumber(stderr,v1); fprintf(stderr,"\n");
              fprintf(stderr,"AValue2 and ABS");  WriteNumber(stderr,a2); WriteNumber(stderr,v2); fprintf(stderr,"\n");
                 * */
            }
            for (j = 0; j < cone.d; j++)
            {
                LinearComb(NewRay[j], Ptr1.RRay[j], v2, Ptr2.RRay[j], v1);
            }
            if (localdebug)
            {
                /*
              fprintf(stderr,"CreatNewRay: New ray ="); WriteArow(stderr, NewRay, cone.d);
                 * */
            }
            Normalize(cone.d, NewRay);
            if (localdebug)
            {
                /*
              fprintf(stderr,"CreatNewRay: Normalized ray ="); WriteArow(stderr, NewRay, cone.d);
                 * */
            }
            AddRay(cone, NewRay);
            //clear(a1); clear(a2); clear(v1); clear(v2);
        }

        public void EvaluateARay1(int i, Cone cone)
        /* Evaluate the ith component of the vector  A x RD.Ray 
            and rearrange the linked list so that
            the infeasible rays with respect to  i  will be
            placed consecutively from First 
         */
        {
            int j;
            double temp, tnext;
            Ray Ptr, PrevPtr, TempPtr;

            temp = 0.0; tnext = 0.0;
            Ptr = cone.FirstRay;
            PrevPtr = cone.ArtificialRay;
            if (PrevPtr.Next != Ptr)
            {
                /*
              fprintf(stderr,"Error.  Artificial Ray does not point to FirstRay!!!\n");
                 * */
            }
            while (Ptr != null)
            {
                temp = 0.0;
                for (j = 0; j < cone.d; j++)
                {
                    tnext = cone.A[i - 1][j] * Ptr.RRay[j];
                    temp += tnext;
                }
                Ptr.ARay = temp;
                /*    if ( temp <= -zero && Ptr != cone.FirstRay) {*/
                if (temp < 0 && Ptr != cone.FirstRay)
                {
                    /* fprintf(stderr,"Moving an infeasible record w.r.t. %ld to FirstRay\n",i); */
                    if (Ptr == cone.LastRay) cone.LastRay = PrevPtr;
                    TempPtr = Ptr;
                    Ptr = Ptr.Next;
                    PrevPtr.Next = Ptr;
                    cone.ArtificialRay.Next = TempPtr;
                    TempPtr.Next = cone.FirstRay;
                    cone.FirstRay = TempPtr;
                }
                else
                {
                    PrevPtr = Ptr;
                    Ptr = Ptr.Next;
                }
            }
            //clear(temp); clear(tnext);
        }

        public void EvaluateARay2(int i, Cone cone)
        /* Evaluate the ith component of the vector  A x RD.Ray 
           and rearrange the linked list so that
           the infeasible rays with respect to  i  will be
           placed consecutively from First. Also for all feasible rays,
           "positive" rays and "zero" rays will be placed consecutively.
         */
        {
            int j;
            double temp, tnext;
            Ray Ptr, NextPtr;
            bool zerofound = false, negfound = false, posfound = false;

            if (cone == null || cone.TotalRayCount <= 0) return;
            temp = 0.0; tnext = 0.0;
            cone.PosHead = null; cone.ZeroHead = null; cone.NegHead = null;
            cone.PosLast = null; cone.ZeroLast = null; cone.NegLast = null;
            Ptr = cone.FirstRay;
            while (Ptr != null)
            {
                NextPtr = Ptr.Next;  /* remember the Next record */
                Ptr.Next = null;     /* then clear the Next pointer */
                temp = 0.0;
                for (j = 0; j < cone.d; j++)
                {
                    tnext = cone.A[i - 1][j] * Ptr.RRay[j];
                    temp += tnext;
                }
                Ptr.ARay = temp;
                /*    if ( temp < -zero) {*/
                if (temp < 0)
                {
                    if (!negfound)
                    {
                        negfound = true;
                        cone.NegHead = Ptr;
                        cone.NegLast = Ptr;
                    }
                    else
                    {
                        Ptr.Next = cone.NegHead;
                        cone.NegHead = Ptr;
                    }
                }
                /*    else if (temp > zero){*/
                else if (temp > 0)
                {
                    if (!posfound)
                    {
                        posfound = true;
                        cone.PosHead = Ptr;
                        cone.PosLast = Ptr;
                    }
                    else
                    {
                        Ptr.Next = cone.PosHead;
                        cone.PosHead = Ptr;
                    }
                }
                else
                {
                    if (!zerofound)
                    {
                        zerofound = true;
                        cone.ZeroHead = Ptr;
                        cone.ZeroLast = Ptr;
                    }
                    else
                    {
                        Ptr.Next = cone.ZeroHead;
                        cone.ZeroHead = Ptr;
                    }
                }
                Ptr = NextPtr;
            }
            /* joining three neg, pos and zero lists */
            if (negfound)
            {                 /* -list nonempty */
                cone.FirstRay = cone.NegHead;
                if (posfound)
                {               /* -list & +list nonempty */
                    cone.NegLast.Next = cone.PosHead;
                    if (zerofound)
                    {            /* -list, +list, 0list all nonempty */
                        cone.PosLast.Next = cone.ZeroHead;
                        cone.LastRay = cone.ZeroLast;
                    }
                    else
                    {                      /* -list, +list nonempty but  0list empty */
                        cone.LastRay = cone.PosLast;
                    }
                }
                else
                {                        /* -list nonempty & +list empty */
                    if (zerofound)
                    {            /* -list,0list nonempty & +list empty */
                        cone.NegLast.Next = cone.ZeroHead;
                        cone.LastRay = cone.ZeroLast;
                    }
                    else
                    {                      /* -list nonempty & +list,0list empty */
                        cone.LastRay = cone.NegLast;
                    }
                }
            }
            else if (posfound)
            {            /* -list empty & +list nonempty */
                cone.FirstRay = cone.PosHead;
                if (zerofound)
                {              /* -list empty & +list,0list nonempty */
                    cone.PosLast.Next = cone.ZeroHead;
                    cone.LastRay = cone.ZeroLast;
                }
                else
                {                        /* -list,0list empty & +list nonempty */
                    cone.LastRay = cone.PosLast;
                }
            }
            else
            {                          /* -list,+list empty & 0list nonempty */
                cone.FirstRay = cone.ZeroHead;
                cone.LastRay = cone.ZeroLast;
            }
            cone.ArtificialRay.Next = cone.FirstRay;
            cone.LastRay.Next = null;
            //clear(temp); clear(tnext);
        }

        public void DeleteNegativeRays(Cone cone)
        /* Eliminate the infeasible rays with respect to  i  which
           are supposed to be consecutive from the head of the Ray list,
           and sort the zero list assumed to be consecutive at the
           end of the list.
         */
        {
            int fii, fiitest;
            double temp;
            Ray Ptr, PrevPtr, NextPtr, ZeroPtr1, ZeroPtr0;
            bool found, completed, zerofound = false, negfound = false, posfound = false;
            bool localdebug = false;

            temp = 0.0;
            cone.PosHead = null; cone.ZeroHead = null; cone.NegHead = null;
            cone.PosLast = null; cone.ZeroLast = null; cone.NegLast = null;

            /* Delete the infeasible rays  */
            PrevPtr = cone.ArtificialRay;
            Ptr = cone.FirstRay;
            if (PrevPtr.Next != Ptr)
            {
                /*
              fprintf(stderr,"Error at DeleteNegativeRays: ArtificialRay does not point the FirstRay.\n");
                 * */
            }
            completed = false;
            while (Ptr != null && !completed)
            {
                /*    if ( (Ptr.ARay) < -zero ){ */
                if (Ptr.ARay < 0)
                {
                    Eliminate(cone, PrevPtr);
                    Ptr = PrevPtr.Next;
                }
                else
                {
                    completed = true;
                }
            }

            /* Sort the zero rays */
            Ptr = cone.FirstRay;
            cone.ZeroRayCount = 0;
            while (Ptr != null)
            {
                NextPtr = Ptr.Next;  /* remember the Next record */
                temp = Ptr.ARay;
                if (localdebug)
                {
                    /*
                    fprintf(stderr,"Ptr.ARay :"); WriteNumber(stderr, temp);
                     * */
                }
                /*    if ( temp < -zero) {*/
                if (temp < 0)
                {
                    if (!negfound)
                    {
                        /*
                      fprintf(stderr,"Error: An infeasible ray found after their removal\n");
                         * */
                        negfound = true;
                    }
                }
                /*    else if (temp > zero){*/
                else if (temp > 0)
                {
                    if (!posfound)
                    {
                        posfound = true;
                        cone.PosHead = Ptr;
                        cone.PosLast = Ptr;
                    }
                    else
                    {
                        cone.PosLast = Ptr;
                    }
                }
                else
                {
                    (cone.ZeroRayCount)++;
                    if (!zerofound)
                    {
                        zerofound = true;
                        cone.ZeroHead = Ptr;
                        cone.ZeroLast = Ptr;
                        cone.ZeroLast.Next = null;
                    }
                    else
                    {/* Find a right position to store the record sorted w.r.t. FirstInfeasIndex */
                        fii = Ptr.FirstInfeasIndex;
                        found = false;
                        ZeroPtr1 = null;
                        for (ZeroPtr0 = cone.ZeroHead; !found && ZeroPtr0 != null; ZeroPtr0 = ZeroPtr0.Next)
                        {
                            fiitest = ZeroPtr0.FirstInfeasIndex;
                            if (fiitest >= fii)
                            {
                                found = true;
                            }
                            else ZeroPtr1 = ZeroPtr0;
                        }
                        /* fprintf(stderr,"insert position found \n %d  index %ld\n",found, fiitest); */
                        if (!found)
                        {           /* the new record must be stored at the end of list */
                            cone.ZeroLast.Next = Ptr;
                            cone.ZeroLast = Ptr;
                            cone.ZeroLast.Next = null;
                        }
                        else
                        {
                            if (ZeroPtr1 == null)
                            { /* store the new one at the head, and update the head ptr */
                                /* fprintf(stderr,"Insert at the head\n"); */
                                Ptr.Next = cone.ZeroHead;
                                cone.ZeroHead = Ptr;
                            }
                            else
                            {                /* store the new one inbetween ZeroPtr1 and 0 */
                                /* fprintf(stderr,"Insert inbetween\n");  */
                                Ptr.Next = ZeroPtr1.Next;
                                ZeroPtr1.Next = Ptr;
                            }
                        }
                        /*
                        Ptr.Next=cone.ZeroHead;
                        cone.ZeroHead=Ptr;
                        */
                    }
                }
                Ptr = NextPtr;
            }
            /* joining the pos and zero lists */
            if (posfound)
            {            /* -list empty & +list nonempty */
                cone.FirstRay = cone.PosHead;
                if (zerofound)
                {              /* +list,0list nonempty */
                    cone.PosLast.Next = cone.ZeroHead;
                    cone.LastRay = cone.ZeroLast;
                }
                else
                {                        /* 0list empty & +list nonempty */
                    cone.LastRay = cone.PosLast;
                }
            }
            else
            {                          /* +list empty & 0list nonempty */
                cone.FirstRay = cone.ZeroHead;
                cone.LastRay = cone.ZeroLast;
            }
            cone.ArtificialRay.Next = cone.FirstRay;
            cone.LastRay.Next = null;
            //clear(temp);
        }

        public void FeasibilityIndices(out int fnum, out int infnum, int i, Cone cone)
        {
            /*Evaluate the number of feasible rays and infeasible rays*/
            /*  w.r.t the hyperplane  i*/
            int j;
            double temp, tnext;
            Ray Ptr;

            temp = 0.0; tnext = 0.0;
            fnum = 0;
            infnum = 0;
            Ptr = cone.FirstRay;
            while (Ptr != null)
            {
                temp = 0.0;
                for (j = 0; j < cone.d; j++)
                {
                    tnext = cone.A[i - 1][j] * Ptr.RRay[j];
                    temp += tnext;
                }
                if (temp >= 0)
                    fnum++;
                else
                    infnum++;
                Ptr = Ptr.Next;
            }
            //clear(temp); clear(tnext);
        }

        public bool LexSmaller(double[] v1, double[] v2, int dmax)
        { /* dmax is the size of vectors v1,v2 */
            bool determined, smaller;
            int j;

            smaller = false;
            determined = false;
            j = 1;
            do
            {
                if (v1[j - 1] != v2[j - 1])
                {  /* 086 */
                    if (v1[j - 1] < v2[j - 1])
                    {  /*086 */
                        smaller = true;
                    }
                    determined = true;
                }
                else
                    j++;
            } while (!(determined) && (j <= dmax));
            return smaller;
        }


        public bool LexLarger(double[] v1, double[] v2, int dmax)
        {
            return LexSmaller(v2, v1, dmax);
        }

        public bool LexEqual(double[] v1, double[] v2, int dmax)
        { /* dmax is the size of vectors v1,v2 */
            bool determined, equal;
            int j;

            equal = true;
            determined = false;
            j = 1;
            do
            {
                if (v1[j - 1] != v2[j - 1])
                {  /* 093c */
                    equal = false;
                    determined = true;
                }
                else
                {
                    j++;
                }
            } while (!(determined) && (j <= dmax));
            return equal;
        }

        public void AddNewHalfspace1(Cone cone, int hnew)
        /* This procedure 1 must be used with PreorderedRun=false 
           This procedure is the most elementary implementation of
           DD and can be used with any type of ordering, including
           dynamic ordering of rows, e.g. MaxCutoff, MinCutoff.
           The memory requirement is minimum because it does not
           store any adjacency among the rays.
        */
        {
            Ray RayPtr0, RayPtr1, RayPtr2, RayPtr2s, RayPtr3;
            int pos1, pos2;
            double prevprogress, progress;
            double value1, value2;
            bool adj, equal, completed;

            value1 = 0.0; value2 = 0.0;
            EvaluateARay1(hnew, cone);
            /*Check feasibility of rays w.r.t. hnew 
              and put all infeasible ones consecutively */

            RayPtr0 = cone.ArtificialRay;   /*Pointer pointing RayPrt1*/
            RayPtr1 = cone.FirstRay;        /*1st hnew-infeasible ray to scan and compare with feasible rays*/
            value1 = cone.FirstRay.ARay;
            if (value1 >= 0)
            {
                if (cone.RayCount == cone.WeaklyFeasibleRayCount) cone.CompStatus = CompStatusType.AllFound;
                return;        /* Sicne there is no hnew-infeasible ray and nothing to do */
            }
            else
            {
                RayPtr2s = RayPtr1.Next;/* RayPtr2s must point the first feasible ray */
                pos2 = 1;
                while (RayPtr2s != null && RayPtr2s.ARay < 0)
                {
                    RayPtr2s = RayPtr2s.Next;
                    pos2++;
                }
            }
            if (RayPtr2s == null)
            {
                cone.FirstRay = null;
                cone.ArtificialRay.Next = cone.FirstRay;
                cone.RayCount = 0;
                cone.CompStatus = CompStatusType.AllFound;
                return;   /* All rays are infeasible, and the computation must stop */
            }
            RayPtr2 = RayPtr2s;   /*2nd feasible ray to scan and compare with 1st*/
            RayPtr3 = cone.LastRay;    /*Last feasible for scanning*/
            prevprogress = -10.0;
            pos1 = 1;
            completed = false;
            while ((RayPtr1 != RayPtr2s) && !completed)
            {
                value1 = RayPtr1.ARay;
                value2 = RayPtr2.ARay;
                CheckEquality(cone.d, RayPtr1, RayPtr2, out equal);
                if (value1 > 0 && value2 < 0 || value1 < 0 && value2 > 0)
                {
                    CheckAdjacency(cone, RayPtr1, RayPtr2, out adj);
                    if (adj) CreateNewRay(cone, RayPtr1, RayPtr2, hnew);
                }
                if (RayPtr2 != RayPtr3)
                {
                    RayPtr2 = RayPtr2.Next;
                    continue;
                }
                if (value1 < 0 || equal)
                {
                    Eliminate(cone, RayPtr0);
                    RayPtr1 = RayPtr0.Next;
                    RayPtr2 = RayPtr2s;
                }
                else
                {
                    completed = true;
                }
                pos1++;
                progress = 100.0 * ((double)pos1 / pos2) * (2.0 * pos2 - pos1) / pos2;
                if (progress - prevprogress >= 10 && pos1 % 10 == 0 && debug)
                {
                    /*
                  fprintf(stderr,"*Progress of iteration %5ld(/%ld):   %4ld/%4ld => %4.1f%% done\n",
                     cone.Iteration, cone.m, pos1, pos2, progress);
                     * */
                    prevprogress = progress;
                }
            }
            if (cone.RayCount == cone.WeaklyFeasibleRayCount) cone.CompStatus = CompStatusType.AllFound;

            //clear(value1); clear(value2);
        }

        public void AddNewHalfspace2(Cone cone, int hnew)
        /* This procedure must be used under PreOrderedRun mode */
        {
            Ray RayPtr0, RayPtr1, RayPtr2;
            Adjacency EdgePtr, EdgePtr0;
            int pos1;
            int fii1, fii2;
            bool localdebug = false;

            EvaluateARay2(hnew, cone);
            /* Check feasibility of rays w.r.t. hnew 
               and sort them. ( -rays, +rays, 0rays)*/

            if (cone.PosHead == null && cone.ZeroHead == null)
            {
                cone.FirstRay = null;
                cone.ArtificialRay.Next = cone.FirstRay;
                cone.RayCount = 0;
                cone.CompStatus = CompStatusType.AllFound;
                return;   /* All rays are infeasible, and the computation must stop */
            }

            if (localdebug)
            {
                pos1 = 0;
                /*
              fprintf(stderr,"(pos, FirstInfeasIndex, A Ray)=\n");
                 * */
                for (RayPtr0 = cone.FirstRay; RayPtr0 != null; RayPtr0 = RayPtr0.Next)
                {
                    pos1++;
                    /*
                  fprintf(stderr,"(%ld,%ld,",pos1,RayPtr0.FirstInfeasIndex);
                     * */
                    /*
                  WriteNumber(stderr,RayPtr0.ARay); 
                  fprintf(stderr,") ");
                     * */
                }
                /*
              fprintf(stderr,"\n");
                 * */
            }

            if (cone.ZeroHead == null) cone.ZeroHead = cone.LastRay;

            EdgePtr = cone.Edges[cone.Iteration];
            while (EdgePtr != null)
            {
                RayPtr1 = EdgePtr.Ray1;
                RayPtr2 = EdgePtr.Ray2;
                fii1 = RayPtr1.FirstInfeasIndex;
                CreateNewRay(cone, RayPtr1, RayPtr2, hnew);
                fii2 = cone.LastRay.FirstInfeasIndex;
                if (fii1 != fii2)
                    ConditionalAddEdge(cone, RayPtr1, cone.LastRay, cone.PosHead);
                EdgePtr0 = EdgePtr;
                EdgePtr = EdgePtr.Next;

                cone.EdgeCount--;
            }
            cone.Edges[cone.Iteration] = null;

            DeleteNegativeRays(cone);

            cone.AddedHalfspaces.Set(hnew, true);

            if (cone.Iteration < cone.m)
            {
                if (cone.ZeroHead != null && cone.ZeroHead != cone.LastRay)
                {
                    if (cone.ZeroRayCount > 200 && debug)
                    {
                        /*
                        fprintf(stderr,"*New edges being scanned...\n");
                         * */
                    }
                    UpdateEdges(cone, cone.ZeroHead, cone.LastRay);
                }
            }

            if (cone.RayCount == cone.WeaklyFeasibleRayCount) cone.CompStatus = CompStatusType.AllFound;

        }


        public void SelectNextHalfspace0(Cone cone, BitArray excluded, out int hnext)
        {
            /*A natural way to choose the next hyperplane.  Simply the largest index*/
            int i;
            bool determined;

            i = cone.m;
            determined = false;
            do
            {
                if (excluded.Get(i))
                    i--;
                else
                    determined = true;
            } while (!determined && i >= 1);
            if (determined)
                hnext = i;
            else
                hnext = 0;
        }

        public void SelectNextHalfspace1(Cone cone, BitArray excluded, out int hnext)
        {
            /*Natural way to choose the next hyperplane.  Simply the least index*/
            int i;
            bool determined;

            i = 1;
            determined = false;
            do
            {
                if (excluded.Get(i))
                    i++;
                else
                    determined = true;
            } while (!determined && i <= cone.m);
            if (determined)
                hnext = i;
            else
                hnext = 0;
        }

        public void SelectNextHalfspace2(Cone cone, BitArray excluded, out int hnext)
        {
            /*Choose the next hyperplane with maximum infeasibility*/
            int i, fea, inf, infmin, fi = 0;   /*feasibility and infeasibility numbers*/

            hnext = 0;

            infmin = cone.RayCount + 1;
            for (i = 1; i <= cone.m; i++)
            {
                if (!excluded.Get(i))
                {
                    FeasibilityIndices(out fea, out inf, i, cone);
                    if (inf < infmin)
                    {
                        infmin = inf;
                        fi = fea;
                        hnext = i;
                    }
                }
            }
            if (debug)
            {
                /*
              fprintf(stderr,"*infeasible rays (min) =%5ld, #feas rays =%5ld\n", infmin, fi);
                 * */
            }
        }

        public void SelectNextHalfspace3(Cone cone, BitArray excluded, out int hnext)
        {
            /*Choose the next hyperplane with maximum infeasibility*/
            int i, fea, inf, infmax, fi = 0;   /*feasibility and infeasibility numbers*/
            bool localdebug = debug;

            hnext = 0;

            infmax = -1;
            for (i = 1; i <= cone.m; i++)
            {
                if (!excluded.Get(i))
                {
                    FeasibilityIndices(out fea, out inf, i, cone);
                    if (inf > infmax)
                    {
                        infmax = inf;
                        fi = fea;
                        hnext = i;
                    }
                }
            }
            if (localdebug)
            {
                /*
              fprintf(stderr,"*infeasible rays (max) =%5ld, #feas rays =%5ld\n", infmax, fi);
                 * */
            }
        }

        public void SelectNextHalfspace4(Cone cone, BitArray excluded, out int hnext)
        {
            /*Choose the next hyperplane with the most unbalanced cut*/
            int i, fea, inf, max, tmax, fi = 0, infi = 0;
            /*feasibility and infeasibility numbers*/

            hnext = 0;

            hnext = 0;

            max = -1;
            for (i = 1; i <= cone.m; i++)
            {
                if (!excluded.Get(i))
                {
                    FeasibilityIndices(out fea, out inf, i, cone);
                    if (fea <= inf)
                        tmax = inf;
                    else
                        tmax = fea;
                    if (tmax > max)
                    {
                        max = tmax;
                        fi = fea;
                        infi = inf;
                        hnext = i;
                    }
                }
            }
            if (!debug)
                return;
            if (max == fi)
            {
                /*
              fprintf(stderr,"*infeasible rays (min) =%5ld, #feas rays =%5ld\n", infi, fi);
                 * */
            }
            else
            {
                /*
              fprintf(stderr,"*infeasible rays (max) =%5ld, #feas rays =%5ld\n", infi, fi);
                 * */
            }
        }

        public void SelectNextHalfspace5(Cone cone, BitArray excluded, out int hnext)
        {
            /*Choose the next hyperplane which is lexico-min*/
            int i, minindex;
            double[] v1, v2;

            minindex = 0;
            v1 = null;
            for (i = 1; i <= cone.m; i++)
            {
                if (!excluded.Get(i))
                {
                    v2 = cone.A[i - 1];
                    if (minindex == 0)
                    {
                        minindex = i;
                        v1 = v2;
                    }
                    else if (LexSmaller(v2, v1, cone.d))
                    {
                        minindex = i;
                        v1 = v2;
                    }
                }
            }
            hnext = minindex;
        }


        public void SelectNextHalfspace6(Cone cone, BitArray excluded, out int hnext)
        {
            /*Choose the next hyperplane which is lexico-max*/
            int i, maxindex;
            double[] v1, v2;

            maxindex = 0;
            v1 = null;
            for (i = 1; i <= cone.m; i++)
            {
                if (!excluded.Get(i))
                {
                    v2 = cone.A[i - 1];
                    if (maxindex == 0)
                    {
                        maxindex = i;
                        v1 = v2;
                    }
                    else if (LexLarger(v2, v1, cone.d))
                    {
                        maxindex = i;
                        v1 = v2;
                    }
                }
            }
            hnext = maxindex;
        }

        public void UniqueRows(int[] OV, int p, int r, double[][] A, int dmax, BitArray preferred, out int uniqrows)
        {
            /* Select a subset of rows of A (with range [p, q] up to dimension dmax) by removing duplicates.
               When a row subset preferred is nonempty, those row indices in the set have priority.  If
               two priority rows define the same row vector, one is chosen.
               For a selected unique row i, OV[i] returns a new position of the unique row i. 
               For other nonuniqu row i, OV[i] returns a negative of the original row j dominating i.
               Thus the original contents of OV[p..r] will be rewritten.  Other components remain the same.
               *uniqrows returns the number of unique rows.
           */
            int i, iuniq, j;
            double[] x;
            bool localdebug = false;

            uniqrows = 0;

            if (p <= 0 || p > r) return;
            iuniq = p; j = 1;  /* the first unique row candidate */
            x = A[p - 1];
            OV[p] = j;  /* tentative row index of the candidate */
            for (i = p + 1; i <= r; i++)
            {
                if (!LexEqual(x, A[i - 1], dmax))
                {
                    /* a new row vector found. */
                    iuniq = i;
                    j = j + 1;
                    OV[i] = j;    /* Tentatively register the row i.  */
                    x = A[i - 1];
                }
                else
                {
                    /* rows are equal */
                    if (preferred.Get(i) && !preferred.Get(iuniq))
                    {
                        OV[iuniq] = -i;  /* the row iuniq is dominated by the row i */
                        iuniq = i;  /* the row i is preferred.  Change the candidate. */
                        OV[i] = j;  /* the row i is tentatively registered. */
                        x = A[i - 1];
                    }
                    else
                    {
                        OV[i] = -iuniq;  /* the row iuniq is dominated by the row i */
                    }
                }
            }
            uniqrows = j;
            if (localdebug)
            {
                /*
              printf("The number of unique rows are %ld\n with order vector",*uniqrows);
                 * */
                for (i = p; i <= r; i++)
                {
                    /*
                  printf(" %ld:%ld ",i,OV[i]);
                     * */
                }
                /*
              printf("\n");
                 * */
            }
        }

        public int Partition(int[] OV, int p, int r, double[][] A, int dmax)
        {
            double[] x;
            int i, j, ovi;

            x = A[OV[p] - 1];

            i = p - 1;
            j = r + 1;
            while (true)
            {
                do
                {
                    j--;
                } while (LexLarger(A[OV[j] - 1], x, dmax));
                do
                {
                    i++;
                } while (LexSmaller(A[OV[i] - 1], x, dmax));
                if (i < j)
                {
                    ovi = OV[i];
                    OV[i] = OV[j];
                    OV[j] = ovi;
                }
                else
                {
                    return j;
                }
            }
        }

        public void QuickSort(int[] OV, int p, int r, double[][] A, int dmax)
        {
            int q;

            if (p < r)
            {
                q = Partition(OV, p, r, A, dmax);
                QuickSort(OV, p, q, A, dmax);
                QuickSort(OV, q + 1, r, A, dmax);
            }
        }

        /*
#ifndef RAND_MAX 
#define RAND_MAX 32767 
#endif
         * */

        public void RandomPermutation(int[] OV, int t, int seed)
        {
            int k, j, ovj;
            double u, xk, r, rand_max = 32767; //(double) RAND_MAX;
            bool localdebug = false;

            Random rnd = new Random(seed);  //******

            for (j = t; j > 1; j--)
            {
                /*
              r=rand();
              u=r/rand_max;
                 * */
                u = rnd.NextDouble();
                xk = (double)(j * u + 1);
                k = (int)xk;
                if (localdebug)
                {
                    /*
                    fprintf(stderr, "u=%g, k=%ld, r=%g, randmax= %g\n", u, k, r, rand_max);
                     * */
                }
                ovj = OV[j];
                OV[j] = OV[k];
                OV[k] = ovj;
                if (localdebug)
                {
                    /*
                    fprintf(stderr, "row %ld is exchanged with %ld\n", j, k);
                     * */
                }
            }
        }

        public void ComputeRowOrderVector(Cone cone)
        {
            int i, itemp;

            cone.OrderVector[0] = 0;
            switch (cone.HalfspaceOrder)
            {
                case RowOrderType.MaxIndex:
                    for (i = 1; i <= cone.m; i++) cone.OrderVector[i] = cone.m - i + 1;
                    break;

                case RowOrderType.MinIndex:
                    for (i = 1; i <= cone.m; i++) cone.OrderVector[i] = i;
                    break;

                case RowOrderType.LexMin:
                case RowOrderType.MinCutoff:
                case RowOrderType.MixCutoff:
                case RowOrderType.MaxCutoff:
                    for (i = 1; i <= cone.m; i++) cone.OrderVector[i] = i;
                    RandomPermutation(cone.OrderVector, cone.m, cone.rseed);
                    QuickSort(cone.OrderVector, 1, cone.m, cone.A, cone.d);
                    break;

                case RowOrderType.LexMax:
                    for (i = 1; i <= cone.m; i++) cone.OrderVector[i] = i;
                    RandomPermutation(cone.OrderVector, cone.m, cone.rseed);
                    QuickSort(cone.OrderVector, 1, cone.m, cone.A, cone.d);
                    for (i = 1; i <= cone.m / 2; i++)
                    {   /* just reverse the order */
                        itemp = cone.OrderVector[i];
                        cone.OrderVector[i] = cone.OrderVector[cone.m - i + 1];
                        cone.OrderVector[cone.m - i + 1] = itemp;
                    }
                    break;

                case RowOrderType.RandomRow:
                    for (i = 1; i <= cone.m; i++) cone.OrderVector[i] = i;
                    RandomPermutation(cone.OrderVector, cone.m, cone.rseed);
                    break;

            }
        }

        public void UpdateRowOrderVector(Cone cone, BitArray PriorityRows)
        /* Update the RowOrder vector to shift selected rows
        in highest order.
        */
        {
            int i, j, k, j1 = 0, oj = 0;
            int rr;
            bool found, localdebug = false;

            if (debug) localdebug = true;
            found = true;
            rr = PriorityRows.Count;
            if (localdebug)
            {
                /*
                set_fwrite(stderr, PriorityRows);
                 * */
            }
            for (i = 1; i <= rr; i++)
            {
                found = false;
                for (j = i; j <= cone.m && !found; j++)
                {
                    oj = cone.OrderVector[j];
                    if (PriorityRows.Get(oj))
                    {
                        found = true;
                        if (localdebug)
                        {
                            /*
                            fprintf(stderr, "%ldth in sorted list (row %ld) is in PriorityRows\n", j, oj);
                             * */
                        }
                        j1 = j;
                    }
                }
                if (found)
                {
                    if (j1 > i)
                    {
                        /* shift everything lower: ov[i].cone.ov[i+1]..ov[j1-1].cone.ov[j1] */
                        for (k = j1; k >= i; k--) cone.OrderVector[k] = cone.OrderVector[k - 1];
                        cone.OrderVector[i] = oj;
                        if (localdebug)
                        {
                            /*
                          fprintf(stderr,"OrderVector updated to:\n");
                             * */
                            for (j = 1; j <= cone.m; j++)
                            {
                                /*
                                fprintf(stderr, " %2ld", cone.OrderVector[j]);
                                 * */
                            }
                            /*
                          fprintf(stderr,"\n");
                             * */
                        }
                    }
                }
                else
                {
                    /*
                  fprintf(stderr,"UpdateRowOrder: Error.\n");
                     * */
                    return;
                }
            }
        }

        public void SelectPreorderedNext(Cone cone, BitArray excluded, out int hh)
        {
            int i, k;

            hh = 0;
            for (i = 1; i <= cone.m && hh == 0; i++)
            {
                k = cone.OrderVector[i];
                if (!excluded.Get(k)) hh = k;
            }
        }

        public void SelectNextHalfspace(Cone cone, BitArray excluded, out int hh)
        {
            if (cone.PreOrderedRun)
            {
                if (debug)
                {
                    /*
                  fprintf(stderr,"debug SelectNextHalfspace: Use PreorderNext\n");
                     * */
                }
                SelectPreorderedNext(cone, excluded, out hh);
            }
            else
            {
                if (debug)
                {
                    /*
                  fprintf(stderr,"debug SelectNextHalfspace: Use DynamicOrderedNext\n");
                     * */
                }

                switch (cone.HalfspaceOrder)
                {

                    case RowOrderType.MaxIndex:
                        SelectNextHalfspace0(cone, excluded, out hh);
                        break;

                    case RowOrderType.MinIndex:
                        SelectNextHalfspace1(cone, excluded, out hh);
                        break;

                    case RowOrderType.MinCutoff:
                        SelectNextHalfspace2(cone, excluded, out hh);
                        break;

                    case RowOrderType.MaxCutoff:
                        SelectNextHalfspace3(cone, excluded, out hh);
                        break;

                    case RowOrderType.MixCutoff:
                        SelectNextHalfspace4(cone, excluded, out hh);
                        break;

                    default:
                        SelectNextHalfspace0(cone, excluded, out hh);
                        break;
                }
            }
        }

#if nemkell
bool Nonnegative(double val)
{
/*  if (val>=-zero) return true;  */
  if (cmp(val,minuszero)>=0) return true;
  else return false;
}

bool Nonpositive(double val)
{
/*  if (val<=zero) return true;  */
  if (cmp(val,zero)<=0) return true;
  else return false;
}

bool Positive(double val)
{
  return !Nonpositive(val);
}

bool Negative(double val)
{
  return !Nonnegative(val);
}

bool EqualToZero(double val)
{
  return (Nonnegative(val) && Nonpositive(val));
}

bool Nonzero(double val)
{
  return (Positive(val) || Negative(val));
}

bool Equal(double val1,double val2)
{
  return (!Larger(val1,val2) && !Smaller(val1,val2));
}

bool Larger(double val1,double val2)
{
  double temp;
  bool answer;

  init(temp);
  sub(temp,val1, val2);
  answer=Positive(temp);
  clear(temp);
  return answer;
}

bool Smaller(double val1,double val2)
{
  return Larger(val2,val1);
}

void abs(double absval, double val)
{
  if (Negative(val)) neg(absval,val);
  else set(absval,val); 
}
#endif

        public void LinearComb(double lc, double v1, double c1, double v2, double c2)
        /*  lc := v1 * c1 + v2 * c2   */
        {
            double temp;

            temp = 0.0;
            lc = v1 * c1;
            temp = v2 * c2;
            lc += temp;
        }

        public void InnerProduct(out double prod, int d, double[] v1, double[] v2)
        {
            double temp;
            int j;
            bool localdebug = false;

            temp = 0.0;
            prod = 0.0;
            for (j = 0; j < d; j++)
            {
                temp = v1[j] * v2[j];
                prod += temp;
            }
            if (localdebug)
            {
                /*
              fprintf(stderr,"InnerProduct:\n"); 
              WriteArow(stderr, v1, d);
              WriteArow(stderr, v2, d);
              fprintf(stderr,"prod ="); 
              WriteNumber(stderr, prod);
              fprintf(stderr,"\n");
                 * */
            }
        }

        /* end of cddcore.c */




    }
}
